
program define rdtab, eclass
	version 14
	syntax varlist(min=3 max=4) [if] [, est_command(string) OLDmatrix(string) NEWmatrix(string) Tau(real 0.0) Alpha(real 0.05) Power(real 0.8) onesided init_val(real 0.001)]
	local y_var = word("`varlist'", 1)
	local treat_var = word("`varlist'", 2)
	local rv_var = word("`varlist'", 3)
	local clust_var = word("`varlist'", 4)
	matrix b = J(1,7,-9999)

	if "`est_command'"=="reg" {
		matrix b[1,1] = _b[treat]
		matrix b[1,2] = _se[treat]
		matrix b[1,3] = e(N)
		if `tau'!=0 matrix b[1,6] = 1 - normal(`tau'/_se[treat] + invnormal(1-(`alpha'/2))) + normal(`tau'/_se[treat] - invnormal(1-(`alpha'/2)))
		if `tau'!=0 & "`onesided'"!="" matrix b[1,6] = normal(`tau'/_se[treat] - invnormal(1-`alpha'))
		qui distinct `clust_var' if e(sample)
		matrix b[1,4] = r(ndistinct)
		gen __xx__absrv = abs(`rv_var')
		qui sum __xx__absrv if e(sample)
		matrix b[1,5] = r(max)
		drop __xx__*
	}

	if "`est_command'"=="rdrobust" {
		matrix b[1,1] = e(tau_cl)
		matrix b[1,2] = e(se_tau_cl)
		matrix b[1,5] = e(h_l)
		matrix b[1,3] = e(N_h_l) + e(N_h_r)
		if `tau'!=0 matrix b[1,6] = 1 - normal(`tau'/e(se_tau_cl) + invnormal(1-(`alpha'/2))) + normal(`tau'/e(se_tau_cl) - invnormal(1-(`alpha'/2)))
		if `tau'!=0 & "`onesided'"!="" matrix b[1,6] = normal(`tau'/e(se_tau_cl) - invnormal(1-`alpha'))
		preserve
		cap keep `if'
		qui distinct `clust_var' if abs(`rv_var') <= e(h_l) & `y_var'!=.
		matrix b[1,4] = r(ndistinct)
		restore
	}
	
	scalar se = b[1,2]
	scalar z = invnormal(1 - `alpha'/2)
	if "`onesided'"!="" scalar z = invnormal(1 - `alpha')
	scalar power = `power'
	scalar init_val = `init_val'
	
	mata: mde_opt = optimize_init()
	mata: optimize_init_which(mde_opt, "max")
	if "`onesided'"=="" mata: optimize_init_evaluator(mde_opt, &compute_hypoth_power())
	if "`onesided'"!="" mata: optimize_init_evaluator(mde_opt, &compute_hypoth_power())
	else mata: optimize_init_evaluator(mde_opt, &compute_hypoth_power_one_sided())
	mata: optimize_init_params(mde_opt, st_numscalar("init_val"))
	mata: optimize_init_argument(mde_opt, 1, st_numscalar("se"))
	mata: optimize_init_argument(mde_opt, 2, st_numscalar("z"))
	mata: optimize_init_argument(mde_opt, 3, st_numscalar("power"))
	mata: mde = optimize(mde_opt)
	mata: st_local("mde", strofreal(mde))
	
	matrix b[1,7] = `mde'
	
	if ("`oldmatrix'"=="" & "`newmatrix'"=="") {
		noi di "Need a new or old matrix name"
	}
	else if ("`oldmatrix'"!="" & "`newmatrix'"!="") {
		noi di "Cannot have old matrix name and new matrix name"
	}
	else if "`oldmatrix'"!="" {
		matrix `oldmatrix' = (`oldmatrix' , b')
	}
	else {
		matrix `newmatrix' = b'
	}
	ereturn post b
end

mata: 
	void compute_hypoth_power(todo, t, s, z, p, y, g, H) {
		y = -(1 - normal(t/s + z) + normal(t/s - z) - p)^2
	}
end

mata:
	void compute_hypoth_power_one_sided(todo, t, s, z, p, y, g, H) {
		y = -(normal(t/s - z) - p)^2
	}
end

